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The Heisenberg ferromagnet (uniaxially anisotropic along z-direction), in the presence of 
time dependent (but uniform over space) magnetic field, is studied by Monte Carlo simulation. 
The time dependent magnetic field was taken as elliptically polarised in such a way that the 
resulting field vector rotates in the XZ-plane. In the limit of low anisotropy, the dynamical 
responses of the system are studied as functions of temperature and the amplitudes of the 
magnetic field. As the temperature decreases, it was found that the system undergoes multiple 
dynamical phase transitions. In this limit, the multiple transitions were studied in details and 
the phase diagram for this observed multicritical behavior was drawn in field amplitude and 
temperature plane. The natures (continuous/discontinuous) of the transitions are determined 
by the temperature variations of fourth-order Binder cumulant ratio and the distributions 
of order parameters near the transition points. The transitions are supported by finite size 
study. The temperature variations of the variances of dynamic order parameter components 
(for different system sizes) indicate the existence of diverging length scale near the dynamic 
transition points. The frequency dependences of the transition temperatures of the multiple 
dynamic transition are also studied briefly. 

INTRODUCTION 

The nonequilibrium dynamical phase transition [1], particularly in the kinetic Ising model, 
has drawn immense interest of the researchers working in the field of nonequilibrium statistical 
physics. The dynamic transition in the kinetic Ising model in the presence of sinusoidally 
oscillating magnetic field was noticed [2] first in the mean field solution of the dynamical 
equation for the average magnetisation. Due to the absence of fluctuations, in the mean field 
study, the dynamic transition would be possible also even in the static (zero frequency) limit. 
But the true dynamic transition should never occur in the static limit due to the presence 
of nontrivial thermodynamic fluctuations. The occurence of the true dynamic transitions for 
models, incorporating the thermodynamic fluctuations, was later shown in several Monte Carlo 
studies [1]. 
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Since, the Ising model is quite simple and a prototype to study the nonequilibrium dynam- 
ical phase transitions, a considerable amount of work was done on it [1]. Among these, several 
studies were done to establish that all these transitions have some features which are similar 
to well known thermodynamic phase transitions. The divergence of "time scale" [3] and the 
"dynamic specific heat" [3] and the divergence of "length scale" [4] are two very important ob- 
servations to establish that the dynamic transition, in kinetic Ising model driven by sinusoidally 
oscillating magnetic field, is indeed a thermodynamic phase transition. 

Although, the dynamic phase transition in kinetic Ising model is an interesting phenomenon 
which acts as a simple example to grasp the various features of nonequilibrium phase transitions, 
it has several limitations. Since the orientations of the spins are limited by only two (up/down) 
directions, some interesting features of dynamic transitions (related to the dynamical transverse 
ordering) cannot be observed in Ising model. The classical vector spin model [5] would be the 
proper choice to study such interesting phenomena which are missing in Ising model. The 
"off-axial" dynamic transition [6], recently observed in Heisenberg ferromagnet, is an example 
of such interesting phenomena. In the "off-axial" dynamic transition, the dynamical symmetry 
along the axis of anisotropy can be broken dynamically by applying an oscillating magnetic 
field along any perpendicular direction. There are several important studies done on classical 
vector spin models driven by oscillating magnetic field. The dynamic phase transition in an 
anisotropic XY ferromagnet driven by oscillating magnetic field was recently studied [7] by 
solving Ginzburg-Landau equation. The dynamic phase transition and the dependence of its 
behaviour on the bilinear exchange anisotropy of classical Heisenberg ferromagnet (thin film), 
was studied by [8] Monte Carlo Simulation. The dynamic transitions alongwith hysteresis 
scaling in Heisenberg ferromagnet was also studied [9] both by Monte Carlo simulation and 
mean field solution. 

All these studies (discussed in above paragraph) of dynamic transitions in classical vector 
spin model give single dynamic transitions. There are a few studies in Heisenberg ferromagnets 
where the multiple dynamic transition was reported [10, 11]. A recent study [10] of dynami- 
cal phase transition in thin Heisenberg ferromagnetic films with bilinear exchange anisotropy 
has shown multiple phase transitions for the surface and bulk layers of the film at different 
temperatures. However, another study [11] shows triple dynamic transitions (at three different 
temperatures), for bulk dynamic ordering only, in uniaxially anisotropic Heisenberg ferromagnet 
driven by elliptically polarised magnetic field. In the present paper, the multiple dynamic tran- 
sition (bulk only) [11] in uniaxially anisotropic Heisenberg ferromagnets driven by elliptically 
polarised magnetic field, is studied extensively. 

In this paper, the multiple dynamic transition, in the uniaxially anisotropic classical Heisen- 
berg ferromagnet in presence of elliptically polarised magnetic field, is studied by Monte Carlo 
simulation. The various dynamical phases were found. The multiple dynamic phase transi- 
tion was observed by studying the temperature variation of the dynamic specific heat and the 
derivatives of dynamic order parameter components. The nature (continuous/discontinuous) 
of the transition was detected by studying the distribution of order parameter and the tem- 
perature variation of Binder cumulant ratio near the transition point. A finite size study was 
done supporting the multiple dynamic phase transitions. A length scale was found to diverge 
near the dynamic transition points. And finally and most importantly, the phase diagram for 
this multiple dynamic phase transition was drawn. The frequency dependence of these multiple 
dynamic transitions was also studied briefly. 
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The paper is organized as follows: in the next section the description of the system, i.e., the 
Heisenberg ferromagnet in polarised magnetic field, is given, section III describes the Monte 
Carlo simulation technique employed here to study the dynamical steady state of this type 
of classical vector spin model, in detail. The numerical results with diagrams are reported in 
section IV and the paper ends with a summary in section V. 



HEISENBERG FERROMAGNET IN POLARIZED MAGNETIC FIELD 

The classical anisotropic (uniaxial and single-site) Heisenberg model, with nearest neighbour 
ferromagnetic interaction in the presence of a magnetic field can be described by the following 
Hamiltonian: 

H = - JE <0> Si • Sj - D^(S iz ) 2 - h ■ S,Si. (1) 

In the above expression, Si[Si X , Si y , Si Z ] represents a classical spin vector (situated at the ith 
lattice site) of magnitude unity, i.e., |S| = 1 or Sf x + Sf y + Sf z = 1. The classical spin vector Si 
may take any (unrestricted) angular orientation in the vector spin space. The first term, in the 
Hamiltonian, represents the nearest -neighbour (< ij >) ferromagnetic (J > 0) interaction. The 
factor D > 0, in the second term, represents the strength of uniaxial (z axis here) anisotropy 
which is favoring the spin vector to be aligned along the z axis. This is readily seen since 
the second term minimizes energy for maximization of the value of Sf. Here, it may be noted 
that for D — > oo the system goes to the Ising limit and for D = the system is in the 
isotropic Heisenberg limit. The last term stands for the interaction with the externally applied 
time dependent magnetic field [h(h x ,h y ,h z )]. The magnetic-field components are sinusoidally 
oscillating in time, i.e., h a = /i 0a cos(c<jt), where h 0a is the amplitude and u is the angular 
frequency (ou = 27r/; / is the frequency) of the ath component of the magnetic field. In the 
present study, the field is taken elliptically polarised. In general, the field vector is represented 
as 

h = xh x + yh y + zh z = xh 0x cos(ut) + zh 0z sm(cut). (2) 
The time eliminated relation between h x (— h 0x cos(ut)) and h z {— /i 02 sin(u;t) is 

w + w = 1 - ® 

This is the equation of an ellipse which indicates that the magnetic field vector is elliptically 
polarised (for ho x ^ h 0z ) and lies on the X — Z plane. For h§ x = h^ z = h (say), the field will 
be circularly polarised and h x + h 2 z = h^. The magnetic fields and the strength of anisotropy 
D are measured in units of J. The model is defined on a simple cubic lattice of linear size L 
with periodic boundary conditions applied in all three (x, y, z) directions. 



MONTE CARLO SIMULATION METHOD 

Monte Carlo (MC) simulation method was employed to study the above described model. 
The algorithm is [12] described below. The system is slowly cooled down from a random initial 
spin configuration [13], to obtain the steady state spin configuration at a particular temperature 



3 



T (measured in units of J/ks, where ks is the Boltzmann constant). The initial random spin 
configuration was generated as follows [6, 13]: two different (uncorrelated) random numbers r 1 
and r 2 (uniformly distributed between -1 and +1), are chosen in such a way that R 2 = r\ + r\ 
becomes less than or equal to unity. The set of values of r\ and r 2 , for which R 2 > 1, are rejected. 
Now, — 2uri, S iy = 2ur 2 and S iz — 1 — 2i? 2 , where u = y/l — R 2 . After preparing initial 
random configuration of spins (this is the proper choice of the spin configurations, corresponding 
to very high temperature), one has to find a steady state configuration for any fixed temperature 
T. For any fixed set of values of ho x , ho z , w and D and at any particular temperature T, a lattice 
site % has been chosen randomly (random updateing scheme). The value (random direction) of 
the spin vector at this randomly chosen site is Si (say). The energy of the system is given by 
the Hamiltonian (Eq. 1). A test spin vector Sj is then chosen (for a trial move) at any random 
direction (following the same algorithm described above). For this choice of Sj, at site % the 
energy will be H = — JS < jj > S / i • Sj — DHi(S' iz ) 2 — h ■ EjSj. The change in energy, associated 
with this change in the direction of spin vector (from Si to Sj at lattice site i), is AH = H' — H. 
The Monte Carlo method [12] will decide whether this trial move is acceptable or not. The 
probability of this move (chosen here) is given by Metropolis rate [14] 

A M 

W(S,-S' I ) = M3n[l,exp(-— )]. (4) 

Now, this probability will be compared with a random number R p (say) (uniformly distributed 
between zero and one). If Rp does not exceed W, the move (Si — > Sj) will be accepted. In this 
way the spin vector Si is updated. L 3 numbers of such random updates of spins, defines one 
Monte Carlo step per site (MCSS) and this is considered as the unit of time in this simulation. 
The efficiency of the MC technique can be increased by choosing the direction of the spin for 
the trial move close to the present direction which will increase the probability of acceptance. 
So, 50 MCSS are required to have one complete cycle of the oscillating magnetic field (of 
frequency / = 0.02) and hence the time period (r) of the field becomes 50 MCSS. Any time 
dependent macroscopic quantity (i.e., any component of magnetisation at any instant t) is 
calculated as follows: Starting with an initial random spin configuration (corresponding to 
the high-temperature disordered phase), the system is allowed to become stable (dynamically) 
up to 5 x 10 4 MCSS (i.e., 1000 complete cycle of the oscillating field). The average value 
of various physical quantities are calculated from further 5 x 10 4 MCSS (i.e., averaged over 
another 1000 cycles). This was checked carefully that the number of MCSS mentioned above is 
sufficient to achieve dynamical steady state value of the measurable quantities, etc. which would 
clearly show the dynamic transition points within limited accuracy. But to describe the critical 
behaviour very precisely (i.e., to estimate critical exponents etc) a much longer run is required. 
The total length of the simulation becomes 10 5 MCSS. The system is slowly cooled down (T has 
been reduced by a small interval 5T = 0.02 here) to get the values of the statistical quantities 
in the lower temperature phase. Here, the last spin configuration corresponding to the previous 
temperature is employed to act as initial configuration for the new (lower) temperature. The 
CPU time required for 10 5 MCSS is approximately 30 min on an Intel Pentium-Ill processor. 

One important point may be noted here regarding the dynamics chosen in this simulational 
study. Since the spin components do not commute with the Heisenberg Hamiltonian, it has 
intrinsic quantum mechanical dynamics. Considering the intrinsic dynamics, there was a study 
[15] of the structure factor and transport properties in XY model. However, the present pa- 
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per aims to study the nonequilibrium phase transition governed by thermal fluctuations only. 
Keeping this in mind, one should choose the dynamics that arising solely only due to the inter- 
action with thermal bath. Since the objective is different, in this paper, the dynamics chosen 
(arise solely due to the interaction with a thermal bath) were Metropolis dynamics. The effect 
of intrinsic spin dynamics is therefore not considered here. 

NUMERICAL RESULTS 

The Monte Carlo simulations, in the present study, were done on a simple cubic lattice 
of linear size L = 20. For a fixed set of values of amplitudes (h 0x ,h 0z ) and frequency (/) of 
polarised magnetic field, strength of anisotropy (D) and temperature (T), the instantaneous 
magnetisation components (per lattice site) were calculated as follows: m x (t) = EjyS'f/L 3 , 
m y (t) = Y^iSf/L 3 and m z (t) = S^f/L 3 . The dynamic order parameter is defined as the time 
averaged magnetisation over a full cycle of the oscillating magnetic field. In this case the 
dynamic order parameter Q is a vector (since the magnetisation is vector). The components 
of the dynamic order parameter are calculated as Q x — (1/t) §m x (t)dt, Q y = (l/T)§m y (t)dt 
and Q z = (1/t) § m z (t)dt. Here, four different dynamical phases are identified. The high 
temperature disordered phase is denoted as Po : (Qx = Q y = Qz = 0). Other three ordered 
phases are P 1 : (Q x = 0, Q y ^ 0,Q Z = 0), P 2 : (Q x ^ 0, Q y = 0,Q Z ^ 0) and P 3 : (Qx = 
0, Q y = 0, Q z 7^ 0). The phase diagram, obtained from multiple dynamic transitions, is plotted 
in h 0x — T plane and is shown in fig 1. The methods of getting the phase boundaries are 
discussed in the following paragraphs. 

The signature of such successive phase transitions are also observed by studying the tem- 
perature variations of the magnitude of the dynamic order parameter, energy, specific heat 
and the temperature derivatives of dynamic order parameter components. For the same set 
of values D = 0.2, h 0x = 0.3 and h 0z = 1.0 the temperature variations of various dynamic 
quantities are studied and shown in Fig. 2. The temperature variations of the dynamic order 
parameter components Q x , Q y and Q z are shown in Fig. 2(a). As the system is cooled down, 
from a high-temperature dynamically disordered (Q = 0) phase, it was observed that first the 
system undergoes a transition from dynamically disordered Po phase (Q = 0) to a dynamically 
Y-ordered (Q y ^ only) phase P 1 and the transition temperature is T cl (say). It may be noted 
here that the resultant vector of elliptically polarised magnetic field lies in the x — z plane and 
the dynamic ordering occurs along the y-direction only (Q x = 0, Q y ^ 0, Q z = 0). This is 
clearly an off-axial transition [6]. In the case of this type of off-axial transition the dynamical 
symmetry (in any direction; y here) is broken by the application of the magnetic field in the 
perpendicular (lies in the x — z plane here) direction. As the system cools down further, this 
phase Pi persists over a considerable range of temperature and at a temperature T c2 a second 
transition was observed. In this phase, the system becomes dynamically ordered both in the 
X and Z directions at the cost of Y ordering. Usually one gets the dynamically ordered phase 
of second kind P 2 : (Q x ^ 0, Q y = 0,Q Z ^ 0). Here, the dynamical ordering is planar (lies on 
x — z plane) and the dynamical ordering occurs in the same plane on which the field vector 
lies. This transition is not of off-axial type. As the temperature decreases further one ends up 
with a low-temperature dynamically ordered phase of third kind P3 : (Q x = 0, Q y = 0, Q z 7^ 0) 
via a transition occurs at temperature T c3 . The system continues to increase the dynamical 
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Z ordering (Q z ^ only) as the temperature decreases further. The low-temperature phase 
is only dynamically Z ordered. No further transition was observed as one cools the system 
further. From Fig. 2 (a), it is quite clear, if one observes carefully, that the sizes of the errorbars 
of Q x , Q y and Q z are maximal near the transition points, indicating the growth of fluctuations 
near the transiton points. 

One gets qualitative ideas of multiple dynamic phase transitions from the above mentioned 
studies. However, to estimate precisely the transition temperatures T ci , T c2 and T c3 further 
studies are required. For this reason, the temperature variations of the derivatives (with respect 
to temperature) of the dynamic order parameter components are studied. The derivatives were 
calculated numerically by using central difference formula [16] 

df(x) = f{x + Sx) - f(x - 5x) 
dx 2Sx U 

Fig. 2(b) shows such variations studied as a function of temperature. The derivative dQ y /dT, 
shows a sharp minimum at T = T cl ~ 1.22. The second transition temperature T c2 was esti- 
mated from the position of sharp maximum of dQ y /dT and the corresponding sharp minima 
of dQ x /dT and dQ z /dT. This gives T c2 ~ 0.96. At slightly lower temperature than T c2 , one 
observes another maximum of dQ x /dT and minimum of dQ z /dT both at the same position 
T = T c3 ~ 0.88. From this study one gets the quantitative measure of the transition temper- 
atures of multiple dynamic transition [11]. The maximum error, in estimating the transition 
temperatures, is 0.01. 

Similar values of the transition temperatures (T cl , T c2 , T c3 ) for the multiple dynamic transi- 
tion can be estimated independently from the study of the temperature variations of dynamic 
energy and specific heat. The dynamical energy (E) is defined as the time averaged value of 
the instantaneous energy over a full cycle of the oscillating magnetic field. From the definition, 
E — (1/t) § Hdt (H is given in Eqn. 1). The temperature variation of E is studied and shown 
in Fig. 2(c). This shows two inflection points and one (the middle one) discontinuity at the 
same location of transition temperatures estimated and mentioned above. This will be very 
clear if one studies the derivative of the dynamic energy, namely the dynamic specific heat 
C = dE/dT (calculated by using central difference formula 5). The temperature variation 
of dynamic specific heat C was studied and shown in Fig. 2(d). It indicates three peaks at 
three different temperatures supporting T ci ~ 1.22, T c2 ~ 0.96 and T c3 ~ 0.88. Thus the 
estimation of transition temperatures for the multiple dynamic transition was reexamined by 
another independent study. The study of the temperature variation of dynamic specific heat 
has another importance. It independently supports multiple transitions as well as it indicates 
that these transitions are indeed thermodynamic phase transitions. Now one may employ both 
the methods of studying the temperature variations of derivatives of dynamic order parameter 
components and the dynamic specific heat to estimate the transition points. One gets three 
different dynamic transitions for the parameter values D = 0.2, h 0x = 0.3 and h 0z = 1.0. This 
three transitions senario is observed for a range of values of h$ x (keeping other parameters fixed 
D = 0.2 and ho z = 1.0) between ho x = 0.1 to ho x = 0.5 (within the precision considered here). 
The temperature variations of dynamic specific heat C are shown in fig 3 for two different 
values of ho x (= 0.1 and 0.5). Both show the three dynamic transitions. It may be noted from 
fig. 3 that the transition temperature decreases as the amplitude h 0x increases. Now, let us see 
what happens if one takes the value of the x-amplitude of elliptically polarised magnetic field, 
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i.e., ho x , outside this range, keeping all other parameter values unchanged. 

For h 0x = 0.0 (effectively the field is now linearly polarised along z direction only), to 
estimate the transition points the temperature variations of dynamic specific heat was studied 
and plotted in fig. 4(a). This shows two distinct and well separated peaks indicating two 
dynamic transitions, one at T ~ 1.24 and other at T ~ 0.98. Now to characterise the different 
phases the temperature variations of dynamic order parameter components were studied. This 
study is shown in fig. 4(b). This shows that the system starts to get dynamically ordered 
from a high-temperature disordered phase. So, the high-temperature ordered phase (Q x ^ 0, 
Qy 7^ 0, Q z = 0) is quite different from Pi (described above for Hq x ^ 0). The second (low- 
temperature ordered phase) is P 3 type. This transition (from high-temperature ordered phase 
to low-temperature ordered phase) occurs at T c3 ~ 0.98. So, one gets two distinct transitions for 
D = 0.2, h 0x = 0.0 and h 0z = 1.0. It may be noted here that the temperature variations of the 
order parameter components, mainly Q x and Q y , are quite scattered (near the low-temperature 
transition). Here, Q x and Q y are calculated by averaging over 2000 cycles discarding first 2000 
cycles (10 5 MCSS). For this reason the derivatives dQ x /dT and dQ y /dT do not give smooth 
variation with respect to temperature and are not shown. 

Now if the value of ho x is higher (say around ho x = 0.6) with all other parameter values fixed 
(i.e., / = 0.02, D = 0.2 and ho z = 1.0) the triple dynamic transitions (observed for ho x = 0.3) 
reduces to double transitions. It gives two distinct peaks of C indicating two different dynamic 
transitions at T ~ 1.20 and T ~ 0.82. So, in this case one gets two dynamically ordered phases, 
namely P 1 (higher temperature ordered phase) and P 3 (lower temperature ordered phase). The 
transition from Pi phase to P 3 phase occurs at T ~ 0.82. The high-temperature ordered phase 
(Pi) transition (from dynamically disordered phase) occurs at T ~ 1.20. These two transitions 
were reexamined by studying the temperature variations of the derivatives of the dynamic order 
parameter and obtained the same transition temperatures. 

From the above discussion it is quite clear that if one studies the dynamic phase transitions 
by varying h 0x only (keeping D = 0.2 and ho z = 1-0 fixed), one would get two transitions 
for hox = 0.0. Just above h§ x = 0.0, i.e., starting from ho x = 0.1 upto ho x = 0.5 one would 
get three transitions (and three ordered phases). Above this value, say h§ x = 0.6 one would 
get again two transitions (two ordered phases). If one continues further, it was observed that 
the two transitions feature continues. The transition points (temperatures) are estimated by 
studying specific heat and the derivatives of the dynamic order parameter components. These 
are not shown in figure, only the results (obtained from the peak positions of the specific heat 
and positions of maximum or minimum of dQ a /dT) are taken. It was observed further that 
all the transition points shift towards lower temperatures for higher values of h 0x . The whole 
results of these multiple transitions temperatures and plotted in T — h 0x plane (fig. 1) to get 
the multiple dynamic phase boundary. 

In the phase diagram (fig. 1), the outermost boundary separtes the dynamically V-ordered 
phase Pi : (Q x = 0, Q y ^ 0, Q z = 0) from the disordered phase Pq : (Q x = 0, Q y = 0, Q z — 0). 
The transition temperature decreases as the value of the field amplitude h§ x increases. For 
lower values of h 0x , the region bounded by the boundaries marked by the symbols circle and 
bullet is the ordered phase characterised as P 2 : (Q x ^ 0, Q y = 0, Q z ^ 0). This region is quite 
small (in area) but very clear and was observed very carefully. The transition temperatures for 
both the boundaries (left and right) of this phase decrease as the field amplitude h§ x increases. 
The innermost (in the T — h 0x plane) region, whose boundary is partly marked by circle (for 
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higher values of ho x ) and partly marked by bullet (for lower values of ho x ), represents the low- 
temperature phase characterised by P 3 : (Q x = 0,Q y = 0,Q Z ^ 0). Here also, the transition 
temperature decreases as the value of the field amplitude h§ x increases. 

The question which should arise now, what will be the nature (continuous/discontinuous) 
of these multiple dynamic transitions ? To get the answer of this questions, the distributions 
of the magnitudes of the dynamic order parameter components are studied near the transition 
points. On the phase boundary (Fig.l), a particular value of h 0x (= 0.7) is chosen. This choice 
is not arbitrary. In this region, the phase boundaries are well separated to study comfortably 
the temperature variations of fourth-order Binder cumulant ratio [17] U y — 1— < Qy > /(3 < 
Qy > 2 ). Since the transition from P to Pi phase is indicated by Q y the distribution of the 
magnitude of Q y is studied near the transition point. The histogram is shown in Fig. 5 for four 
different temperatures around the transition temperature. This figure shows that as one goes 
through the results of temperature from T = 1.16 to that of higher temperature, the peak of 
the distributions shifts towards Q y = continuously and around T = 1.20 this gets peaked 
near Q y ~ 0. For slightly higher value of the temperature, say T = 1.22, the distribution is 
peaked around Q y = 0, corresponding to the disordered phase. This indicates the transition 
is continuous with T cl ~ 1.20 (for the transition from P to P\ phase) and one could guess it 
from the temperature variation of Q y and dQ y /dT. This same technique was applied to study 
the nature of the transition from Pi to P 3 phase. The distribution of Q y is shown (in Fig.6) for 
three different temperatures. Here, the different (from the previous case just discussed above) 
senario was observed and may be noted, that the value of Q y drops to zero from a nonzero value 
as the system gets cooled. At T = 0.80 the distribution is singly peaked near Q y ~ 0.63. At 
slightly lower temperature T = 0.79, the distribution get doubly peaked. One additional peak 
appears near Q y ~ 0. This simultaneous appearance of two peaks at very close to the transition 
temperature is a clear signature of discontinuous transition with T c ~ 0.79 (for the transition 
from Pi to P 3 phase). From this study, one can get the idea about the maximum error (0.01 
here) in estimation of the transition temperature. As one lowers the temperatures, T = 0.78, 
one gets only one peak of the distribution around Q y = after the transition. Here also, this 
could be guessed from the temperature variation of Q y . These studies indicate that the high- 
temperature transition is a continuous one and the low-temperature transition is discontinuous 
type. These natures of the transitions were reexamined by studying the temperature variation 
of Binder cumulant ratio U y . For a continuous transition, U y should change monotonically 
from to 2/3 as one tunes the system from the disordered to the ordered phase. On the other 
hand, for a discontinuous transition, U y develops sharp minimum, whose location corresponds 
to the transition point. The results, obtained here, indeed show this (Fig. 7). Fig. 7 shows, 
that as one tunes the system from high temperature, the Binder cumulant ratio U y first grows 
monotonically from to 2/3 around T cl ~ 1.20 indicating that the high-temperature transition 
is continuous. As the temperature decreases further, it shows a very sharp minimum around 
T c ~ 0.80 (for the transition from Pi to P 3 phase) indicating the discontinuous nature of the 
low-temperature transition. These two (multi-histogram analysis and Binder cumulant ratio) 
studies together confirm that the high-temperature transition is a continuous type and the 
low-temperature transition is discontinuous type. Since, the evidence of the nature of the high- 
temperature transition is found to be of continuous type, it will be quite interesting to study 
the scaling behaviour (if any) with the estimations of critical exponents. But for this, one has 
to estimate the transition temperatures more precisely which will require larger sizes of the 
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system. 

What will be the nature of the transition from P 2 to P 3 phase ? For lower values of h 0x 
(=0.3 chosen here) one would be able to see this transition. From Fig.2(a), this transition 
is indicated by the temperature variation of Q z . Since, this transition corresponds to very 
small change of Q z , it is quite difficult to observe any change of the distribution of Q z near 
the transition temperature. However, the temperature variation of the Binder cumulant ratio 
U z = 1— < Q\ > /(3 < Ql > 2 ) indeed shows (see Fig. 8) smooth variation near the transition 
point and it grows monotonically from to 2/3 which indicates the transition is continuous 
with T c3 ~ 0.88 (for the transition from P 2 to P3 phase). All these natures of the transitions 
are marked in the phase diagram (Fig.l). 

The finite size study is also done to confirm that the observation is not an artifact of limited 
size of the system. This study was done for a particular set (/ = 0.02, D = 0.2, h 0x = 0.3, h 0z = 
1.0) of values of the parameters. Different statistical quantities were studied as functions of 
temperatures for different linear sizes (L = 10, 20, 30 here) of the system. This particular choice 
of the values of D, h 0x and h 0z is meaningful in the sense that the three transitions phenomenon 
was observed clearly for this particular values of parameters and the most important part of 
the phase boundary. Figure 9 shows the temperature variations of dQ a /dT for L = 10, 20 and 
30. From the figure it is clear that both the sharpness and the height (depth) of maximum 
(minimum) indicating the transition points increase as the system size increase. The dynamic 
specific heat is also studied as a function of temperature for different L and same set of other 
parameter values. This is shown in fig. 10. This also indicates that the height of the peaks 
(indicating the transitions) increases as the system size increases. The finite size study, done 
and briefly reported here, at least may indicate that the multiple dynamic transitions mentioned 
above is not an artifact of finite size effect. 

The variances of the dynamic order parameter components i.e., L 3 Vai(Q a ) = L 3 [< Q 2 a > 
— < Qa > 2 ] = L 3 5Q 2 a are studied as function of temperature and for different system sizes. 
The values of other parameters are / = 0.02, D = 0.2, h 0x = 0.3 and h 0z = 1.0. Figure 11(a) 
shows temperature variation of L 3 Var(Q x ) = L 3 [< Q 2 X > — < Q x > 2 ] = L 3 5Q 2 X for L = 10,20 
and 30. It gets sharply peaked at around T = T c2 ~ 0.96 (the transition temperature from Pi 
phase to P 2 phase). The height of the peak increases for larger system sizes. In fig. 11(b) the 
L 3 Vai(Q y ) is plotted against temperature for different system sizes. It shows two peaks. The 
high-temperature peak (very short in height but visible in this scale) is located at T ~ 1.22. 
This corresponds to the transition from disordered (Po) to first Y-ordered (Pi) phase. The low- 
temperature peaks are positioned at around T ~ 0.96 corresponding to the transition from Pi 
to P 2 . Here, also the heights of the peaks increases systematically as the system size increases. 
Similar things are observed for L 3 Vai(Q z ) and shown in fig. 11(c). Here, an important thing 
should be noted that this study indicates that their exists a diverging length scale at the 
dynamic transition points [4], for the multiple dynamic transitions in this system. It should 
also be noted that the transition (from P 2 to P3) temperature T C 3 cannot be resoved from T c2 
in the present study. It was already indicated by the results shown in Fig. 2 (a). The sizes of 
the errorbars of Q x , Q y and Q z becomes maximum near T cl and T c2 . However, the growth of 
the errorbars of Q z is not quite clear near T c3 . 

The frequency dependence of these multiple dynamic transitions are also studied very briefly. 
For a lower frequency (/ = 0.01) of the polarised magnetic field, the multiple dynamic transition 
was studied with D = 0.2, h 0x = 0.3 h 0z = 1.0. The temperature variations of the magnitudes 
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of the dynamic order parameter components are plotted in Fig. 12(a). This shows similar 
kind of temperature variations as that was obtained for / = 0.02 (D = 0.2, h 0x = 0.3 and 
h 0z = 1.0), with considerable shifting of transition points towards lower temperatures. From 
the study of the temperature variations of the derivatives of the dynamic order parameter 
components (shown in Fig.l2(b)) one gets T cl ~ 1.06, T c2 ~ 0.76 and T c3 ~ 0.62. A similar 
study was done for higher values of h 0x (= 1.5) (/ = 0.01, D = 0.2 and h 0z = 1.0). The 
results are shown in Fig. 13. In this case only two transitions were observed and the transition 
temperatures estimated are T c \ ~ 0.60 and T C 3 ~ 0.30. From these studies, one can get the 
indication that as the frequency decreases (or as one approaches the static limit) the transition 
points shifts towards lower temperatures keeping the qualitative nature of the phase diagram 
(for / = 0.02) unchanged. This is quite expected result, since as frequency decreases the 
magnetisation components can follow the field in time and the dynamic transition disappears. 
This is why the transition is called truely dynamic and discussed in details in earlier studies [1] 
in the case of Ising model. 

SUMMARY 

The classical Heisenberg ferromagnet (uniaxially anisotropic) in presence of a time-varying 
polarised magnetic field and in contact with a thermal bath is studied by Monte Carlo simula- 
tion using Metropolis dynamics. The magnetic field is elliptically polarised and the field vector 
rotates on XZ - plane. For a very weakly anisotropic (uniaxial) system and particular set of pa- 
rameter values (amplitudes of field in x and z directions) if the system cools down it undergoes 
successive phase transitions. For a range of values of field amplitude along x-direction, three 
dynamic phase transitions were observed. Keeping the values of anisotropy strength and am- 
plitude of field along z direction fixed, the phase transitions were studied for different values of 
the field amplitude along the x direction. In a plane formed by the temperature and the ratio of 
field amplitudes the phase diagram for the multiple dynamic phase transition was drawn. The 
nature (continuous/discontinuous) of the transitions are investigated by the temperature vari- 
ations of the Binder cumulant ratio and the distribution of order parameter near the transition 
points. The finite size study was done to check that the transitions are not artifacts of limited 
system size. The variances of dynamic order parameter components are studied as a function 
of temperature taking system size as parameter. This particular study indicates the existence 
of a diverging length scale at the transition points. The choice of the parameter values, in the 
present study, are obtained from various trials. In other parameter range the multiple transi- 
tions are also possible to observe, however this particular choice gives quite distinct results of 
the observed multicritical behaviour [18]. Since, the transition across the outermost boundary 
of the phase diagram is found to be continuous, it will be much interesting to study the scaling 
behaviour and to estimate the critical exponents. The frequency variations of these multiple 
dynamic transitions are also studied very briefly. This shows that the qualitative nature of the 
phase diagram remains unchanged however the transition temperatures are lowered by lowering 
the frequency of the polarised magnetic field. This is, of course, an expected result experienced 
from the earlier studies [1], particularly, in the case of the dynamic transition in Ising model. 

The dynamic transition in Ising ferromagnet can be explained simply by spin reversal and 
nucleation [19]. The multiple dynamic transition occurs in Heisenberg ferromagnet possibly 
due to the coherent spin rotation [13]. To get the clear idea about it one has to study also the 
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dynamic configurations of spins in details. The study of the relaxation dynamics of the system 
in the low anisotropic limit may help to analyse few results. 

The alternative methods of studying this multiple dynamic phase transitions in anisotropic 
Heisenberg ferromagnet driven by polarised magnetic field are: (i) to study Landau-Lifshitz- 
Gilbert equation of motion [20] with Langevin dynamics, (ii) spin wave analysis of anisotropic 
Heisenberg ferromagnet in presence of polarised magnetic field. 
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Figure 1: The phase diagram in T — h$ x plane for the multiple dynamic transitions for D = 0.2 
/ = 0.02 and h 0z = 1.0. C.T. stands for continuous transition and D.T. stands for discontinuous 
transition. The natures (continuous/discontinuous) of the transitions across the different phase 
boundary are marked by C.T and D.T. 

Figure 2: The temperature variations of different dynamical quantities for / = 0.02, D = 0.2, 
h-ox = 0.3 and h Qz = 1.0. (A) The magnitudes of components of dynamic order parameter Q x , 
Q y and Q z . Vertical lines are errorbars. Note that the sizes of the errorbars become maximum 
near the transition points. (B) derivatives of dynamic order parameters, (C) the dynamic 
energy and the (D) dynamic specific heat. Solid lines are guide to the eye. The vertical arrows 
in (C) and (D) indicates the transition points. The maximum error in estimating the transition 
temperature is 0.01. 

Figure 3: The temperature variations of dynamic specific heat for / = 0.02, D = 0.2, h 0z = 1.0. 
(A) for h 0x = 0.1 and (B) for h 0x = 0.5. Solid lines are guide to the eye. The maximum error 
in estimating the transition temperature is 0.01. 

Figure 4: The temperature variations of (A) dynamic specific heat and (B) magnitudes of the 
order parameter components. For / = 0.02, D = 0.2, h 0x = 0.0 and h 0z = 1.0. Vertical lines are 
errorbars. Note that the sizes of the errorbars becomes maximum near the transition points. 
Solid line in (A) is guide to the eye. 

Figure 5: The unnormalised distributions of the magnitudes of dynamic order parameter com- 
ponent Q y for different temperatures. Here, / = 0.02, D = 0.2, h 0x = 0.7, h 0z = 1.0 and 
/ = 0.02. 

Figure 6: The unnormalised distributions of the magnitudes of the dynamic order parameter 
component Q y for different temperatures. Here, D = 0.2, ho x = 0.7, h 0z = 1.0 and / = 0.02. 

Figure 7: The temperature variation of the Binder cumulant ratio U y near the transition points. 
Here, D = 0.2, h 0x = 0.7, h 0z = 1.0 and / = 0.02. 

Figure 8: The temperature variation of the Binder cumulant ratio U z near the transition point. 
Here, D = 0.2, h 0x = 0.3, h 0z = 1.0 and / = 0.02. 

Figure 9: Temperature variations of derivatives of dynamic order parameter components for 
different system sizes (L = 10,20 and 30) for / = 0.02, D = 0.2, h 0x = 0.3 and h 0z = 1.0. (A) 
dQ x /dT versus T. (B) dQ y /dT versus T. (C) dQ z /dT versus T. Solid lines are guide to the 
eye. 

Figure 10: Temperature variations of dynamic specific heat for different system sizes [L = 10, 
20 and 30) for / = 0.02, D = 0.2, h 0x = 0.3 and h 0z = 1.0. Solid lines are guide to the eye. 
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Figure 11: Temperature variations of variances of the dynamic order parameter components 
for different system sizes (L = 10, 20 and 30) for D = 0.2, h 0x = 0.3 and h 0z = 1.0. (A) L 3 5Q 2 X 
versus T. (B) L 3 5Q 2 versus T. (C) L 3 5Q 2 Z versus T. Solid lines are guide to the eye. 



Figure 12: Temperature variations of the (A) magnitudes of dynamic order parameter com- 
ponents (B) the derivatives of the order parameter components. Here, D = 0.2, ho x = 0.3, 
h 0z = 1.0 and / = 0.01. 



Figure 13: Temperature variations of the (A) magnitudes of dynamic order parameter com- 
ponents (B) the derivatives of the order parameter components. Here, D = 0.2, h 0x = 1.5, 
h 0z = 1.0 and / = 0.01. 
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